https://gsalvatovallverdu.gitlab.io/python/curve_fit/

Isochoric Nucleation Detection analysis | stage 3: Poisson modelling¶

As part of the PhD work of Bruno M. Guerreiro © 2024. If using this notebook, please cite the paper: https://doi.org/10.1021/acsbiomaterials.2c00075

Disclaimer: due to the changing nature of programming documentation, lab work developed and tacit knowledge in this notebook, please contact the author at bruno.guerreiro@fulbrightmail.org if something is not working properly. The code is not actively maintained.

General considerations¶

  • Temperature data must be positive
  • The unfrozen fraction must be in the natural log form

Preparing the data¶

In [1]:
# manage data and fit
import pandas as pd
import numpy as np
import plotly.graph_objs as go
import matplotlib.pyplot as plt

# first part with least squares
from scipy.optimize import curve_fit

# second part about ODR
from scipy.odr import ODR, Model, Data, RealData

# style and notebook integration of the plots
import seaborn as sns
%matplotlib inline
sns.set(style="whitegrid")
In [2]:
# Import data
%store -r water_Tnuc
%store -r fp025_Tnuc
%store -r fp05_Tnuc
%store -r chi_water
%store -r chi_fp025
%store -r chi_fp05

water = pd.DataFrame()
fp025 = pd.DataFrame()
fp05 = pd.DataFrame()

water['chi'] = np.log(pd.Series(chi_water))
fp025['chi'] = np.log(pd.Series(chi_fp025))
fp05['chi'] = np.log(pd.Series(chi_fp05))

water['temp'] = pd.Series(water_Tnuc) *-1
fp025['temp'] = pd.Series(fp025_Tnuc) *-1
fp05['temp'] = pd.Series(fp05_Tnuc) *-1
In [3]:
### Manually choosing what datapoints to fit


# Water
Twater_manFit = water.temp[120:len(water.temp)-10]
Chiwater_manFit = water.chi[120:len(water.chi)-10]

# Fucopol 0.25%
Tfp025_manFit = fp025.temp[30:len(fp025.temp)-25]
Chifp025_manFit = fp025.chi[30:len(fp025.chi)-25]

# Fucopol 0.5%
Tfp05_manFit = fp05.temp[15:len(fp05.temp)-15]
Chifp05_manFit = fp05.chi[15:len(fp05.chi)-15]

Model fitting¶

In [4]:
import math

def poisson_model(T, gamma, n):
    
    ### CONSTANTS ###
    A = 2.9
    beta = -2
    Tm = 0
    
    # division of constants
    const = A/beta
    # numerator
    num = np.multiply(gamma, np.power(np.subtract(T, Tm), 1+n))
    # denominator
    denom = 1+n
    # do division
    tmp = const * (num/denom)
    # exponentiate
    return tmp
In [5]:
Chi_fit_water = []
for i in Twater_manFit:
    Chi_fit_water.append(poisson_model(i, 0.0001, 4))
    
Chi_fit_fp025 = []
for i in Tfp025_manFit:
    Chi_fit_fp025.append(poisson_model(i, 0.0001, 4))
    
Chi_fit_fp05 = []
for i in Tfp05_manFit:
    Chi_fit_fp05.append(poisson_model(i, 0.0001, 4))
In [6]:
# Plot data
from plotly.subplots import make_subplots
fig = make_subplots(rows=1, cols=3)

# Experimental data
fig.add_trace(go.Scatter(name='Water', x=water.temp, y=water.chi, mode='markers', marker_color='rgba(0,0,0,.2)'), row=1, col=1)
fig.add_trace(go.Scatter(name='FP 0.25%', x=fp025.temp, y=fp025.chi, mode='markers', marker_color='rgba(0,0,255,.2)'), row=1, col=2)
fig.add_trace(go.Scatter(name='FP 0.5%', x=fp05.temp, y=fp05.chi, mode='markers', marker_color='rgba(255,0,0,.2)'), row=1, col=3)

# First-approximation fittings
fig.add_trace(go.Scatter(name='Water', x=Twater_manFit, y=Chi_fit_water, mode='lines', line_color='rgba(0,0,0,1)'), row=1, col=1)
fig.add_trace(go.Scatter(name='FP 0.25%', x=Tfp025_manFit, y=Chi_fit_fp025, mode='lines', line_color='rgba(0,0,255,1)'), row=1, col=2)
fig.add_trace(go.Scatter(name='FP 0.5%', x=Tfp05_manFit, y=Chi_fit_fp05, mode='lines', line_color='rgba(255,0,0,1)'), row=1, col=3)


fig.update_layout(template='simple_white', height=400, showlegend=False,
                 xaxis2_title='<b>-Tnuc (ºC)<b>',
                 yaxis_title='<b>ln(Chi)<b>',
                 font_family='Arial', font_color='black', font_size=16
                 )

fig.show()

Fit optimization¶

In [7]:
RMSE_water = []
RMSE_fp025 = []
RMSE_fp05 = []


for i, j in zip(water.chi, Chi_fit_water):
    RMSE_water.append(np.sqrt(np.mean((i-j)**2)))
    
for i, j in zip(fp025.chi, Chi_fit_fp025):
    RMSE_fp025.append(np.sqrt(np.mean((i-j)**2)))
    
for i, j in zip(fp05.chi, Chi_fit_fp05):
    RMSE_fp05.append(np.sqrt(np.mean((i-j)**2)))
In [8]:
from scipy.optimize import curve_fit
popt_water, pcov_water = curve_fit(f=poisson_model, xdata=Twater_manFit, ydata=Chiwater_manFit, p0=(0.01, 1), sigma=RMSE_water, maxfev=12000)
popt_fp025, pcov_fp025 = curve_fit(f=poisson_model, xdata=Tfp025_manFit, ydata=Chifp025_manFit, p0=(0.000001, 10), sigma=RMSE_fp025, maxfev=12000)
popt_fp05, pcov_fp05 = curve_fit(f=poisson_model, xdata=Tfp05_manFit, ydata=Chifp05_manFit, p0=(0.000001, 10), sigma=RMSE_fp05, maxfev=12000)

Results of fit¶

In [9]:
print('DATA FOR WATER:','\nOptimized value for gamma: ',popt_water[0], '\nOptimized value for n: ',popt_water[1])

DATA FOR WATER: 
Optimized value for gamma:  6.4199715065082686e-12 
Optimized value for n:  10.269109653531459
In [10]:
print('DATA FOR 0.25% FUCOPOL:','\nOptimized value for gamma: ',popt_fp025[0], '\nOptimized value for n: ',popt_fp025[1])

DATA FOR 0.25% FUCOPOL: 
Optimized value for gamma:  2.7893101596601606e-15 
Optimized value for n:  14.199805211490776
In [11]:
print('DATA FOR 0.5% FUCOPOL:','\nOptimized value for gamma: ',popt_fp05[0], '\nOptimized value for n: ',popt_fp05[1])

DATA FOR 0.5% FUCOPOL: 
Optimized value for gamma:  3.7948659545168766e-35 
Optimized value for n:  33.37507578001123
In [12]:
# compute a standard deviation error from pcov
gamma_opt_water, n_opt_water = popt_water
gamma_opt_fp025, n_opt_fp025 = popt_fp025
gamma_opt_fp05, n_opt_fp05 = popt_fp05

# Water
perr_water = np.sqrt(np.diag(pcov_water))
Dgamma_water, Dn_water = perr_water

# FucoPol 0.25%
perr_fp025 = np.sqrt(np.diag(pcov_fp025))
Dgamma_fp025, Dn_fp025 = perr_fp025

# FucoPol 0.5%
perr_fp05 = np.sqrt(np.diag(pcov_fp05))
Dgamma_fp05, Dn_fp05 = perr_fp05

print("ERROR DATA FOR WATER:")
print("Gamma = %6.30f +/- %4.30f" % (gamma_opt_water, Dgamma_water))
print("n = %6.3f +/- %4.3f" % (n_opt_water, Dn_water))
print("")
print("ERROR DATA FOR FUCOPOL 0.25%:")
print("Gamma = %6.30f +/- %4.30f" % (gamma_opt_fp025, Dgamma_fp025))
print("n = %6.3f +/- %4.3f" % (n_opt_fp025, Dn_fp025))
print("")
print("ERROR DATA FOR FUCOPOL 0.5%:")
print("Gamma = %6.40f +/- %4.40f" % (gamma_opt_fp05, Dgamma_fp05))
print("n = %6.3f +/- %4.3f" % (n_opt_fp05, Dn_fp05))

ERROR DATA FOR WATER:
Gamma = 0.000000000006419971506508268586 +/- 0.000000000002543765378348991293
n = 10.269 +/- 0.169

ERROR DATA FOR FUCOPOL 0.25%:
Gamma = 0.000000000000002789310159660161 +/- 0.000000000000000624509598625161
n = 14.200 +/- 0.095

ERROR DATA FOR FUCOPOL 0.5%:
Gamma = 0.0000000000000000000000000000000000379487 +/- 0.0000000000000000000000000000000000404543
n = 33.375 +/- 0.446
In [13]:
R2_water = (np.corrcoef(Twater_manFit, Chiwater_manFit)[0,1])**2
R2_fp025 = (np.corrcoef(Tfp025_manFit, Chifp025_manFit)[0,1])**2
R2_fp05 = (np.corrcoef(Tfp05_manFit, Chifp05_manFit)[0,1])**2


print('R2 of water: ', round(R2_water,3))
print('R2 of FP 0.25%: ', round(R2_fp025,3))
print('R2 of FP 0.5%: ', round(R2_fp05,3))

R2 of water:  0.784
R2 of FP 0.25%:  0.841
R2 of FP 0.5%:  0.918
In [14]:
# Plot data with optimized model
fig = go.Figure()

fig = make_subplots(rows=1, cols=3)

# Water
fig.add_trace(go.Scatter(name='Exp', x=water.temp, y=water.chi, mode='markers', marker_color='rgba(0,0,0,.2)'), row=1, col=1)
fig.add_trace(go.Scatter(name='Fit', x=Twater_manFit, y=poisson_model(Twater_manFit, gamma_opt_water, n_opt_water), mode='lines', line_width=15, line_color='rgba(0,0,0,.5)'), row=1, col=1)
fig.add_trace(go.Scatter(name='Fit', x=water.temp, y=poisson_model(water.temp, gamma_opt_water, n_opt_water), mode='lines', line_color='rgba(0,0,0,.8)'), row=1, col=1)

# FP 0.25%
fig.add_trace(go.Scatter(name='Exp', x=fp025.temp, y=fp025.chi, mode='markers', marker_color='rgba(0,0,255,.2)'), row=1, col=2)
fig.add_trace(go.Scatter(name='Fit', x=Tfp025_manFit, y=poisson_model(Tfp025_manFit, gamma_opt_fp025, n_opt_fp025), mode='lines', line_width=15, line_color='rgba(0,0,255,.5)'), row=1, col=2)
fig.add_trace(go.Scatter(name='Fit', x=fp025.temp, y=poisson_model(fp025.temp, gamma_opt_fp025, n_opt_fp025), mode='lines', line_color='rgba(0,0,255,.8)'), row=1, col=2)

# FP 0.5%
fig.add_trace(go.Scatter(name='Exp', x=fp05.temp, y=fp05.chi, mode='markers', marker_color='rgba(255,0,0,.2)'), row=1, col=3)
fig.add_trace(go.Scatter(name='Fit', x=Tfp05_manFit, y=poisson_model(Tfp05_manFit, gamma_opt_fp05, n_opt_fp05), mode='lines', line_width=15, line_color='rgba(255,0,0,.5)'), row=1, col=3)
fig.add_trace(go.Scatter(name='Fit', x=fp05.temp, y=poisson_model(fp05.temp, gamma_opt_fp05, n_opt_fp05), mode='lines', line_color='rgba(255,0,0,.8)'), row=1, col=3)

fig.update_layout(template='simple_white', showlegend=False, height=400,
                 xaxis1_title='',
                 xaxis3_title='',
                 xaxis2_title='<b>-Tnuc (ºC)<b>',
                 xaxis1_range=[4,15],
                 yaxis_title='<b>ln(Chi)<b>',
                 font_family='Arial', font_color='black', font_size=16
                 )

fig.show()
In [15]:
fig = make_subplots(rows=1, cols=3)

# Water
fig.add_trace(go.Scatter(name='Exp', x=np.exp(water.chi), y=water.temp*-1, mode='markers',marker_color='rgba(0,0,0,.2)', marker_size=8, marker_symbol=100),row=1,col=1)
fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(water.temp, gamma_opt_water, n_opt_water)), y=water.temp*-1, mode='lines', line_color='rgba(0,0,0,1)'),row=1,col=1)

# FP 0.25%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp025.chi), y=fp025.temp*-1, mode='markers',marker_color='rgba(0,0,255,.2)', marker_size=8, marker_symbol=100),row=1,col=2)
fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(Tfp025_manFit, gamma_opt_fp025, n_opt_fp025)), y=Tfp025_manFit*-1, mode='lines', line_color='rgba(0,0,255,1)'),row=1,col=2)

# FP 0.5%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp05.chi), y=fp05.temp*-1, mode='markers',marker_color='rgba(255,0,0,.2)', marker_size=8, marker_symbol=100),row=1,col=3)
fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(fp05.temp, gamma_opt_fp05, n_opt_fp05)), y=fp05.temp*-1, mode='lines', line_color='rgba(255,0,0,1)'),row=1,col=3)


fig.update_layout(template='simple_white', width=1000, height=500, showlegend=False,
                 xaxis2_title='<b>Unfrozen fraction (%)<b>',
                 yaxis_title='<b>Nucleation temperature (ºC)<b>',
                 yaxis1_range=[-15,-5],
                 font_family='Arial', font_color='black', font_size=16,
                 xaxis_range=[-0.05,1.05])

fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2)

fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2, tickprefix="<b>",ticksuffix ="</b>")
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1, tickprefix="<b>",ticksuffix ="</b>")


fig.show()
In [16]:
fig = go.Figure()

# Water
fig.add_trace(go.Scatter(name='Exp', x=np.exp(water.chi), y=water.temp*-1, mode='markers',marker_color='rgba(0,0,0,.8)', marker_size=8, marker_symbol=115))
#fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(water.temp, gamma_opt_water, n_opt_water)), y=water.temp*-1, mode='lines', line_color='rgba(0,0,0,.25)'))

# FP 0.25%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp025.chi), y=fp025.temp*-1, mode='markers',marker_color='rgba(0,0,255,.8)', marker_size=8, marker_symbol=15))
#fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(fp025.temp, gamma_opt_fp025, n_opt_fp025)), y=fp025.temp*-1, mode='lines', line_color='rgba(0,0,255,.25)'))

# FP 0.5%
fig.add_trace(go.Scatter(name='Exp', x=np.exp(fp05.chi), y=fp05.temp*-1, mode='markers',marker_color='rgba(255,0,0,.8)', marker_size=8, marker_symbol=15))
#fig.add_trace(go.Scatter(name='Fit', x=np.exp(poisson_model(fp05.temp, gamma_opt_fp05, n_opt_fp05)), y=fp05.temp*-1, mode='lines', line_color='rgba(255,0,0,.25)'))


fig.update_layout(template='simple_white', width=350, height=500, showlegend=False,
                 xaxis_title='<b>Unfrozen fraction (%)<b>',
                 yaxis_title='<b>Nucleation temperature (ºC)<b>',
                 font_family='Arial', font_color='black', font_size=14,
                 xaxis_range=[-0.05,1.05])

fig.update_traces(marker_size=8, marker_symbol=0)

fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2)
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1)

fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2, tickprefix="<b>",ticksuffix ="</b>")
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1, tickprefix="<b>",ticksuffix ="</b>")


fig.show()
#fig.write_image("images/survivor.svg")

Induction time¶

$$IT = 1/J$$$$J = \gamma\Delta{T}^n$$
In [17]:
J_water = np.multiply(gamma_opt_water, np.power(water.temp, n_opt_water))
J_fp025 = np.multiply(gamma_opt_fp025, np.power(fp025.temp, n_opt_fp025))
J_fp05 = np.multiply(gamma_opt_fp05, np.power(fp05.temp, n_opt_fp05))

# Fitting
T_range = np.linspace(0.01,16, 100)
IT_water_fit = 1/np.multiply(gamma_opt_water, np.power(T_range, n_opt_water))
IT_fp025_fit = 1/np.multiply(gamma_opt_fp025, np.power(T_range, n_opt_fp025))
IT_fp05_fit = 1/np.multiply(gamma_opt_fp05, np.power(T_range, n_opt_fp05))

IT_water = 1/J_water
IT_fp025 = 1/J_fp025
IT_fp05 = 1/J_fp05

#####
#IT_waterOptimal_fit = 1/np.multiply(gamma_opt_water, np.power(T_range, n_opt_water))
#IT_waterRepeat_fit = 1/np.multiply(6.23e-12, np.power(T_range, 7.8))
In [18]:
fig = go.Figure()

# Fits
fig.add_trace(go.Scatter(name='Water', x=T_range*-1, y=IT_water_fit, mode='lines', line_width=3, line_color='rgba(0,0,0,.75)'))
fig.add_trace(go.Scatter(name='FucoPol 0.5%', x=T_range*-1, y=IT_fp05_fit, mode='lines', line_width=3, line_color='rgba(255,0,0,.75)'))
fig.add_trace(go.Scatter(name='FucoPol 0.25%', x=T_range*-1, y=IT_fp025_fit, mode='lines', line_width=3, line_color='rgba(0,0,255,.75)'))


fig.update_layout(template='simple_white', width=600, height=500, showlegend=True,
                 yaxis_title='<b>Induction Time<b>',
                 xaxis_title='<b>Working temperature (ºC)<b>',
                 font_family='Arial', font_color='black', font_size=14)

tickvals = [1/60000, 1/60, 1, 60, 1440, 10080, 40320, 483840, 4838400, 48384000, 483840000, 48384000000]
ticktext = ['Millisecond','Second', 'Minute', 'Hour', 'Day', 'Week', 'Month', 'Year', 'Decade', 'Century', 'Millennium', '∞']

fig.update_xaxes(mirror=True, showline=True, ticks='outside', linewidth=2, tickprefix="<b>",ticksuffix ="</b>",
                range=[-4, -13])
fig.update_yaxes(mirror=True, showline=True, ticks='outside', linewidth=2, dtick=1, type='log', 
                 tickprefix="<b>",ticksuffix ="</b>", tickvals=tickvals, ticktext=ticktext, range=[-0.5,7])



fig.show()